Local mean-field study of capillary condensation in silica aerogels 
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We apply local mean-field (i.e. density functional) theory to a lattice model of a fluid in contact 
with a dilute, disordered gel network. The gel structure is described by a diffusion-limited cluster 
aggregation model. We focus on the influence of porosity on both the hysteretic and the equilibrium 
behavior of the fluid as one varies the chemical potential at low temperature. We show that the 
shape of the hysteresis loop changes from smooth to rectangular as the porosity increases and that 
this change is associated to disorder-induced out-of-equilibrium phase transitions that differ on 
adsorption and on desorption. Our results provide insight in the behavior of 4 He in silica aerogels. 

PACS numbers: 64.60.-i,68.45.Da,75.60.Ej 
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I. INTRODUCTION 

The influence of quenched disorder on phase transi- 
tions and critical phenomena continues to be the focus of 
intensive experimental and theoretical research activity. 
Major effects are expected and actually observed when 
the disorder couples linearly to the order parameter of the 
system, a situation that is realized when a fluid or a fluid 
mixture is confined within a porous glass or is in contact 
with the interconnected strands of a gel. A dilute rigid 
network is a particularly interesting medium from a the- 
oretical perspective because exclusion (i.e., confinement) 
effects do not play a dominant role, so that the random- 
field Ising model (RFIM) may be a useful framework to 
interpret the experimental observations 0. 

A striking example of the influence of a gel network 
on fluid phase behavior is provided by the thermody- 
namic studies of Chan and co-workers on 4 He in silica 
aerogels of varying porosity. Silica aerogels are highly 
porous, fractal-like solids made of a tenuous network of 
SiC>2 strands interconnected at random sites. In a 95% 
porosity aerogel, specific heat and adsorption (vapor- 
pressure isotherm) measurements performed in the 
vicinity of the critical temperature of the pure fluid 
(T C =5.195K) show evidence of a phase separation be- 
tween a low-density "vapor" phase (presumably com- 
posed of 4 He vapor plus a liquid film around the silica 
strands 3]) and a high-density "liquid" phase filling the 
whole pore space. The first-order transition appears to 
terminate at a sharply defined critical point that is only 
31mK below T c , which suggests that the system is in 
a weak random-field regime. The coexistence boundary 
in the presence of aerogel is, however, much narrower 
than in the pure system. Subsequently, similar results 
were obtained with N2 in the same aerogel, using light 
scattering and vapor-pressure isotherms[4|. Since out-of- 
cquilibrium and hysteretic effects due to domain forma- 
tion are generally observed in random-field systems be- 
low the critical temperature of the pure system (as illus- 
trated by the behavior of binary mixtures in aerogels 
and diluted antiferromagnets in a magnetic field 6]), it is 
noteworthy that no hysteresis is present in the measure- 



ments performed in Refs.0,0]. The gel- fluid system has 
thus reached equilibrium within the time scale of the ex- 
periments, which is indeed found much longer than the 
characteristic time of activated dynamics (J] ■ The situ- 
ation changes, however, at lower temperatures and the 
adsorption of 4 He in a 98% porosity aerogel is clearly 
hysteretic at 3.60K and 2.34K0, @' Both adsorption 
and desorption isotherms display a vertical step at well- 
defined pressures, but draining occurs at a lower pressure 
than filling|9(. The shape of the hysteresis loop, more- 
over, depends on the porosity of the aerogel: in an aero- 
gel of 87% porosity, 4 He adsorbs and desorbs gradually at 
2.42K and there is no signature of a "liquid- vapor" phase 
coexistence (see Fig. 4(c) in Ref. 8]). Such hysteretic be- 
havior is reminiscent of capillary condensation in a low- 
porosity solid like Vycor glass where one observes a rapid 
increase of the adsorbed quantity at a pressure below the 
bulk saturated vapor pressure, but no sharp vertical step 
in the adsorption isotherms |lCj. The mechanism for hys- 
teresis in porous substrates has prompted much discus- 
sion in the literature and various explanations have been 
proposed, focusing either on single-pore metastability a 
la van der Waals or on network pore-blocking effect sjllj. 
Since both models seem completely inadequate to de- 
scribe light aerogels, the 4 He experiments raise several 
questions: what is the scenario for the change in the 
shape of the hysteresis loops 12]? Do filling and draining 
obey different mechanisms? What is the true equilibrium 
behavior when hysteresis is present? 

These questions are addressed in the present work 
where we build on our earlier studies of capillary con- 
densation in disordered porous solids [H Q HU, but fo- 
cusing on aerogels. We are mainly interested in the influ- 
ence of the porosity on the hysteretic behavior of sorp- 
tion isotherms. This is partly motivated by theoretical 
studies of the zero-temperature RFIM which predict the 
existence of a disorder-induced out-of-equilibrium phase 
transition in the hysteresis loop^fj- The first experi- 
mental observations of this phase transition have been 
recently reported in the literature and we argue that 
the change in the 4 He adsorption isotherm from sharp to 
smooth can be interpreted within the same framework. 
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We propose a different scenario for draining and we re- 
late the observed behavior to an out-of-equilibrium phase 
transition associated to the depinning of the liquid- vapor 
interface. The approach we develop is similar to that 
used previously for calculating the irreversible behavior 
of spin glasses, random-field ferromagnets and diluted 
antiferromagnets|l8j : it consists in studying the evolu- 
tion of the free-energy surface (more precisely, the grand- 
potential surface f2) as the external driving field (here, 
the pressure P of the external vapor or, equivalently, the 
chemical potential fj,) is changed. f2 is a functional of 
the local fluid density and is treated in the mean-field 
approximation. At low temperatures this multidimen- 
sional free-energy landscape is characterized by a large 
number of local minima in which the system may be 
trapped. The main physical assumption underlying our 
description is that thermally activated processes play a 
negligible role on the time scale of the experiments (in 
other words, the dynamics is similar to that at T = 0). 
The evolution of the system then proceeds either contin- 
uously by the deformation of the local minimum in which 
the system is trapped or when this latter becomes unsta- 
ble by a jump to another minimum (avalanche); the re- 
sponse to the driving field is then discontinuous and irre- 
versible. Since it is impossible to perform such a study for 
a continuous model because of computational limitations 
(especially when finite-size scaling analysis is required 
near phase transitions), we adopt a lattice-gas descrip- 
tion that incorporates at a coarse-grained level the geo- 
metric and energetic disorder of the gel- fluid mixture|l9|. 
Previous work has shown that many of the phenomena 
observed in experiments on fluids in disordered porous 
solids can be re pro duced q ualitatively by such simple lat- 
tice models [T3, [3 il5t l20l . As in other studies of phase 
transitions in aerogel[2lH22|. we model the solid by a frac- 
tal structure obtained by diffusion-limited cluster-cluster 
aggregation. 

The paper is organized as follows. In Sec. II, we intro- 
duce the lattice model and discuss some important fea- 
tures of the aerogel structure. In Sec. Ill, we present the 
local mean-field theory and describe the numerical pro- 
cedure. The results for the adsorption, desorption, and 
equilibrium isotherms in 87% and 95% porosity aerogels 
at T/T c = 0.45 are given in Sec. IV. The final section 
presents a summary and conclusions. 



II. LATTICE MODEL AND AEROGEL 
STRUCTURE 

The lattice model used in this work describes the solid 
as a collection of fixed impurities that exert a random yet 
correlated (by the connectivity of the strands) external 
field on the atoms of the fluid. The Hamiltonian is given 
byHHH 




FIG. 1: Three-dimensional realization of a 98% porosity 
DLCA aerogel on a bec lattice of size L = 100 with periodic 
boundary conditions. 
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where n — 0,1 is the usual fluid occupation variable 
(i = 1...N) and r\i — 1,0 is a quenched random variable 
that characterizes the presence of gel particles on the lat- 
tice (when rji = 0, site i is occupied by the gel); Wff and 
w s f denote, respectively, the fluid-fluid and solid-fluid 
attractive interactions, \x is the fluid chemical potential, 
and the double summations run over all distinct pairs of 
nearest-neighbor (n.n.) sites. One can thus vary the gel 
porosity, p = (l/N)J2i r lij by changing the number of 
solid sites or modify the "wettability" of the solid sur- 
face by changing the ratio y = w s f/wjf. For y = 1/2, 
the model reduces to a site diluted Ising model, as can be 
seen by transforming the fluid occupation variable 7$ to 
an Ising spin variable, Si — 2r^ — l|24j: preferential ad- 
sorption of the liquid phase onto the gel is thus modeled 
by y > 1/2. Random fields are generated in the sys- 
tem when y ^ 1/2, and the fluctuating part of the field 
acting on spin Si is proportional to the number of solid 
particles that occupy the nearest neighbors of site i[24|. 
This a discrete random variable that can take the values 
0, 1, ...c, where c is the coordination number of the lattice, 
and whose probability distribution is strongly porosity- 
dependent. 

Gel configurations (i.e., sets {iji}) are generated by a 
standard on-lattice diffusion-limited cluster-cluster ag- 
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FIG. 2: Log-log plot of the aerogel correlation function g(r) 
for p = 0.99, 0.98, 0.97, 0.95, 0.92 and 0.87 (for top to bottom). 
The solid line is a fit with slope —1.13 that corresponds to the 
fractal regime for p — 0.99. 



FIG. 3: Average cluster size £g (estimated from the location 
of the first minimum of g{r)) as a function of the gel concen- 
tration pg — 1 — P- The dashed line corresponds to the power 
law fit £ G = 0.707p G ' 868 . 



gregation (DLCA) algorithm po) adapted to a body- 
centered cubic (bcc) lattice with periodic boundary con- 
ditions. The choice of the bcc lattice will be motivated 
in Sec. IV. B. The DLCA algorithm mimics the growth 
process of base catalyzed aerogels used in helium exper- 
iments, and it has been shown to reproduce the main 
structural features of these aerogels measured from scat- 
tering experiments p6| . A typical exemple of a 98% 
porosity DLCA aerogel on a lattice of linear size L = 100 
is shown in Fig. 1 (from now on, we take the lattice con- 
stant a as the unit length; the total number of sites in 
the lattice is thus N — 2L 3 ). One can clearly see the 
fractal-like character of the gel network that results from 
the aggregation mechanism. 

More quantitative information on the structural prop- 
erties of the gel can be extracted from the two-point 
correlation function g(r). Each curve shown in Fig. 2 
corresponds to a different porosity p and results from 
an average over several simulations in a box of size 
L = 100 (for p = 0.87,0.92,0.95) or L = 200 (for 
p = 0.97,0.98,0.99). Only the most dilute samples ex- 
hibit a true intermediate regime described by the power- 
law behavior g(r) ~ r — ' 3_d f) revealing the fractal char- 
acter of the intra-cluster density-density correlations (2?|. 
For p — 0.99, one finds df w 1.87, a value that lies 
in the range 1.7 — 1.9 expected for DLCA in three di- 
mensions. Note, however, that the fractal dimension de- 
creases with increasing porosity 0, so that the asymp- 
totic value must be somewhat smaller. Another esti- 
mation of df can be obtained from the average cluster 

size £g by assuming that £g varies as Pq 1 ^ 3 , where 
PG = 1 — P is the gel concentration. As suggested in 
Ref. [2(| , £g can be estimated from the location of the first 
minimum of g(r) (which is hardly visible on the scale of 
Fig. 2). The plot of £g versus pc shown in Fig. 3 yields 



dt « 1.85. Since the correlation length is in the range 
650 - 1300A for a 98% base-catalyzed aerogel[!|29| and 
£,g{pg = 0.02) w 21 in the simulation, one can estimate 
that the lattice constant a corresponds to about 30 — 60 A. 
This is consistent with the coarse-grained picture of a gel 
site representing a Si02 particle with a diameter of about 
30l. 

Another relevant length scale is the size of the largest 
cavity in the aerogel 30]. Fig. 4 shows the distribution 
P(n) of nearest distances n from an empty site to the 
aerogel as a function of porosity (the integer "distance" 
n is defined here as the length of the shortest path on 
the lattice from an empty site to a gel site; n = 1 means 
that the empty site is n.n. of a gel site). The plot is the 
normalized histogram of these distances and J^i P( n ') 
gives the probability of being closer than n to the aero- 
gel. One can see that the size of the largest cavity is 
n w 32 for a 99% aerogel and decreases to n rj 10 for a 
95% aerogel and to n » 5 for a 87% aerogel. It is no- 
ticeable that the distribution changes significantly when 
decreasing the porosity from 95% to 87%. In the former 
case, P(n) has its maximum at n — 2 and there is a 
significant proportion of empty sites that are not in the 
immediate vicinity of the gel. In the latter case, P(n) 
is monotonic and strongly peaked at n = 1, which indi- 
cates that most of the empty sites are very close to the 
gel. We shall see in section IV. A that these differences 
have important implications for the behavior of the fluid 
during adsorption. 

It is clear from Fig. 3 that the minimum system size 
necessary to describe correctly collective effects occuring 
inside the aerogel on long length scales (like a sharp con- 
densation event in the whole pore space) depends strongly 
on the porosity. For instance, a box of linear size L = 100 
is not large enough to represent the whole pore space of 
a 99% aerogel because it does not contain a sufficient 
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Specifically, the results presented in this work are calcu- 
lated with y = 2, a value for which the filling of the 87% 
aerogel at the lower closure point of the hysteresis loop 
is roughly the same as in the experiment at T — 2.42K 
(i.e., at T/T c k 0.45). This corresponds to about 1/3 
of the filling reached at the plateau just before satura- 
tion, as can be seen in Fig. 4(c) of Ref.[j| (see also Fig. 
5(a) below). All calculations are done at this single re- 
duced temperature (recall that kT c /wft — c/4 = 2 in 
the mean-field approximation). 



III. LOCAL MEAN-FIELD THEORY (LMFT) 
AND NUMERICAL PROCEDURE 

For the present model, the LMFT consists in solving 
the self-consistent equations for the thermally averaged 
fluid densities Pi({r]i}) =< >t obtained from mini- 
mizing the mean-field grand-potential functional|l3j 



FIG. 4: Distribution P(n) of nearest distances from an 
empty site to the gel (see text). From right to left: p — 
0.99, 0.98, 0.97, 0.95, 0.92, 0.87 (the dashed lines are guides for 
the eye). 



number of connected fractal aggregates (a problem that 
cannot be cured by merely increasing the number of re- 
alizations). Within the framework of local mean- field 
theory, it is however impossible to simulate much larger 
lattices because the computational time and the memory 
storage requirement become prohibitive (at finite temper- 
ature, the local fluid densities are continuous variables, 
which forbids the use of bits algorithms). In the present 
work, we consider 95% and 87% aerogels for which we can 
investigate the statistics of collective events while work- 
ing with lattices of reasonable size (between L = 25 and 
L = 100). 

Since it has been shown previously that the interface 
between the porous solid and the bulk gas may have a 
dramatic effect on the desorption process [l^, two differ- 
ent setups are considered. In the first one, the system 
is periodically replicated in all directions so that no in- 
terfacial effects are taken into account. In the second 
setup an interface is created by placing a slab of vapor of 
width L\, = 10 (the gas "reservoir") in contact with one 
of the [100] faces of the simulation box. Periodic bound- 
ary conditions are then imposed in the y and z directions 
and reflective boundary conditions in the x direction (ob- 
tained by reflecting the lattice at the boundaries). 

In order to completely specify the model, one must also 
fix the value of the interaction parameter y. As shown 
previously |l3, the shape of the hysteresis loops changes 
with y. Since it is meaningless in a coarse-grained picture 
to compute y from the actual values of the fluid-fluid and 
solid-fluid van der Waals interactions, we have chosen its 
value so as to reproduce approximately the height of the 
hysteresis loop in the 87% aerogel at low temperature. 



^({Pi}) = ^BTyjpi In + (rji - pi) ln(?7i - pi)] 
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The variational procedure SCl/Spi = gives a set of N 
coupled non-linear equations 
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where (3 = 1/ (ksT) and the sum runs over the c nearest 
neighbors of site i. At low temperature, these equations 
may have several solutions, and the grand potential cor- 
responding to solution {pf } is given bv[l4| 

fi« = k B TY,Vi ln(l - J) + w ff ^ p?p° . (4) 

i ^ l <ij> 

By using an iterative method to solve Eqs. (3), one only 
finds solutions that are only local minima of the grand- 
potential surface, i.e., metastable states[Is|. 

For a given realization of the aerogel, the sorption 
isotherms (i.e., the curves pf = (1/N) ^ Pi({Vi}> M)) are 
obtained by increasing or decreasing the chemical po- 
tential in small steps Sp. At each subsequent p, the con- 
verged solution at p — Sp (on adsorption) or at p + 5p (on 
desorption) is used to start the iterations. The isotherms 
are then averaged over a number of gel realizations de- 
pending on the system size. In order to determine the 
equilibrium isotherms, we have also searched for addi- 
tional solutions of Eqs. (3) inside the hysteresis loop by 
performing scanning trajectories, as explained below in 
section IV. C. 

To accelerate the convergence of the numerical pro- 
cedure, the iterations have been updated (i.e., the new 




FIG. 5: Hysteresis loops in the 87% (a) and 95% (b) model FIG. 6: Representative adsorption isotherms in 87% (a) and 
aerogels calculated at T/T c = 0.45. Equilibrium isotherms 95% (b) aerogels. System sizes are L = 50 (a) and L = 100 
are indicated by the dotted lines. (b). Notice the change in the scale of the x axis between the 

two figures. 



IV. RESULTS AND DISCUSSION 



value of pi is substituted into the r.h.s. of Eqs. (3) before 
waiting for the rotation through the indices i to be com- 
plete) and the evolution of each site density pi has been 
monitored in the following way. Fluid sites are divided 
into two categories, active and passive. At the beginning 
of the calculation, all sites are active and the procedure 
stops when all sites are passive. At each iteration only 
active sites are considered and their new density is calcu- 
lated using Eqs. (3). If | p[ n) - | < 1CT 6 , 
where n denotes the nth iteration, the site becomes pas- 
sive; otherwise its remains active and its nearest-neighors 
become (or remain) active (some sites can thus be passive 
during a few iterations and become active again). The 
advantage of this algorithm is that the number of active 
sites may quickly become a small fraction of the total 
number of sites, which of course significantly reduces the 
overall computation time. For instance, the number of 
active sites is only 0(L 2 ) when a planar liquid- vapor in- 
terface propagates through the system. The algorithm is 
useful in dilute aerogels because of the presence of large 
cavities. We have checked that the numerical results are 
not different (to the precision of the calculations) from 
those obtained without updating or with using a more 
standard convergence criterion, like in our previous stud- 
ies of fluid adsorption in random matrices|l3l Il4l Il5| . 
Convergence (to an accurracy of 10 -6 ) typically requires 
between 10 2 to 10 3 iterations, and several CPU hours on 
a 2.4 GHz workstation are needed to calculate a single ad- 
sorption isotherm in a system of linear size L — 100 (us- 
ing a step Sp/wff — 10~ 2 or 10~ 3 to increment the chem- 
ical potential). The search for equilibrium isotherms is 
much more time consuming (see below in Sec. IV. C) and 
has been made possible by using parallel computation on 
a Beowulf cluster of 24 processors. 



The results of the present study are summarized in 
Fig. 5 that shows the adsorption, desorption, and equi- 
librium isotherms calculated for the 87% and 95% model 
aerogels for y — 2 and T/T c = 0.45 (from now on, p 
is in units of Wff, the fluid-fluid interaction parameter). 
These curves result from the detailed numerical analysis 
described in the following and correspond to the thermo- 
dynamic limit. The most striking feature is the change 
in the shape of the hysteresis loop from smooth to rect- 
angular as the porosity increases. This change is similar 
to that observed experimentally by Chan and co-workers 
(see Figs. 4(b) and 4(c) in Ref.jg]). As is explained be- 
low, we predict that the three isotherms (adsorption, des- 
orption, and equilibrium) are discontinuous in the 95% 
aerogel (at least within mean-field theory where ther- 
mal fluctuations are neglected). The underlying physical 
mechanisms are, however, quite different. 

A. Adsorption 

We first examine the adsorption process. Calcula- 
tions performed in the presence and in the absence 
of a gel/reservoir interface yield isotherms that are al- 
most indistinguishable (except for the smallest system 
sizes), confirming the conclusion reached in our previous 
work|l5| that adsorption does not depend on the exis- 
tence of a free surface. This is also in line with the re- 
sults of previous calculations done in single slit-like or 
cylindrical pores |3l|. 

In Figs. 6(a) and 6(b) are shown some adsorption 
isotherms calculated in 87% and 95% porosity samples 
of linear size L = 50 and L = 100, respectively (this cor- 
responds to roughly the same ratio L/£q ~ 10). In both 
cases, we focus on the region where the adsorption is the 
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FIG. 7: Cross-section of a 87% aerogel sample on adsorption 
at n = -5,-4.64,-4.5 and -4.35 (L=100). Gel sites are 
shown in black and fluid density is shown as various degrees of 
grey. As explained in the text, one can actually distinguishes 
only two regions corresponding respectively to pi < 0.05 (light 
grey) and pi > 0.95 (dark grey). 



steepest (see Fig. 5). At lower coverages, the isotherms 
look gradual and smooth in the two systems, but as /i in- 
creases, they consist of little steps of varying sizes. In the 
87% aerogel, the size of these "avalanches" remains small 
all the way up to the slowly increasing plateau that ex- 
tends to saturation. On the other hand, in the 95% aero- 
gel, the size of the jumps tends to increase with \i and in 
most of the samples there is a large final avalanche after 
which the filling is almost complete (there are however 
some rare realizations where the last jump is not much 
larger that the preceeding ones, as illustrated by one of 
the isotherms in Fig. 6(b)). 

In order to better visualize the underlying microscopic 
mechanism, we show in Figs. 7 and 8 some cross-sections 
of typical 87% and 95% porosity samples for different val- 
ues of the chemical potential along the isotherms (here, 
the size of both systems is L = 100). One can see that 
in both aerogels the first stage of the adsorption process 
is the formation of a liquid layer that coats the aero- 
gel strands. Then, as \x increases, the film thickens and 
condensation occurs in the smallest crevices defined by 
neighboring gel strands. In the 87% aerogel, it is diffi- 
cult to discriminate between these two filling processes 
because the available empty space is small, as illustrated 
by Fig. 4. For the same reason, the vapor bubbles that 
remain in the system at p, = —4.5 are isolated, and, as 
they shrink in size, the adsorption continues gradually 
until the solid is completely filled with liquid. This is 
precisely the scenario described in Ref.[8j from the ex- 
perimental observations with 4 He. This is also similar 
to what happens in a low-porosity glass like Vvcor|32T|. 
On the contrary, in the 95% sample, one clearly distin- 
guishes the small capillary condensation events that oc- 
cur in some regions of the aerogels (compare the figures 



FIG. 8: Same as Fig. 7 for a 95% sample at fi, = 
-4.5, -4.16, -4.13 and -4.125 (L=100). This is one of the 
samples shown in Fig. 5.b 



for p — —4.16 and p = —4.13) and the major final event 
that corresponds to the filling of a large void space span- 
ning the whole sample (as p, is increased from —4.13 to 
—4.125). Note that the local "liquid- vapor" interfaces in- 
side the gel are very sharp at this low temperature, which 
explains why one can only distinguishes two different re- 
gions in Figs. 7 and 8. Indeed, as illustrated in Fig. 9 
for the lighter aerogel, it is found that the distribution 
of the pi's is bimodal, with most of the fluid sites having 
a density lower than 0.05 or larger than 0.95 (the distri- 
bution is similar in the 87% porosity aerogel with just a 
little more intermediate densities). 

From these results (see for instance the cross-section 
of the 95% aerogel at p, = —4.13), there is no indication 
that the radius of curvature of the liquid- vapor interface 
is concave and uniform throughout the sample just be- 
fore the major condensation event, as was suggested in 
Ref. 8]. This makes unlikely any interpretation in terms 
of a traditional capillary condensation model based for 
instance on the application of the Kelvin equation|3j,|33j. 
Indeed, condensation in the remaining void space is it- 
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FIG. 9: Histogram of the local fluid densities pi in the 95% 
aerogel sample of Fig. 8 at p — —4.13. 
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FIG. 10: Average adsorption isotherms in 87% (a) and 95% 
(b) aerogels for different system sizes. The number of gel 
realizations ranges from 500 (L = 35) to 20 (L = 100) for the 
87% aerogel, and from 2000 (L = 35) to 300 (L = 100) for 
the 95% aerogel 
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FIG. 11: Scaling plot of the compressibility curves dpf/dfi 
during adsorption in the 95% aerogel with \ = 1-22 (Jl t (L) = 
-4.123,-4.126,-4.128 for L = 50,70,100, respectively; ex- 
trapolation to L — > oo gives fj-t =lirnL^oo71 t (Z/) ~ —4.131). 



self triggered by a condensation event that occurs in one 
(or several) small region of the aerogel and that induces 
further condensation in the system. This collective be- 
havior in the form of avalanches of varying sizes is simi- 
lar to that of the Gaussian RFIM at T — proposed by 
Sethna and co-workers to describe the Barkhausen effect 
in low-T ferromagnetic materials Increasing slowly 
the chemical potential of the adsorbed fluid is equivalent 
to increasing adiabatically the applied field in a ferromag- 
net, and changing the porosity of the aerogel modifies the 
amount of quenched disorder in the system. In the RFIM, 
the amount of disorder is controlled by the width of the 
random-field distribution: only small avalanches are ob- 
served for large disorder, resulting in a smooth hysteresis 
loop, whereas one macroscopic avalanche produces a dis- 
continuity in the magnetization for small disorder |l6ll34j . 
The transition between these two regimes for a certain 
value of the disorder corresponds to a critical out-of- 
equilibrium phase transition at which the distribution of 
avalanches follows a power law. 

The existence of a major condensation event in most 
of the 95% samples of size L = 100 therefore suggests 
that there is one macroscopic (i.e., infinite) avalanche 
in the thermodynamic limit, corresponding to a finite 
jump in the adsorption isotherm. However, this can be 
only confirmed by performing a finite-size scaling analy- 
sis of the isotherms obtained after averaging over many 
gel realizations. We indeed observe significant sample-to- 
sample fluctuations both in the location and the height 
of the largest jump. Such average isotherms are shown 
in Fig. 10. For the 87% aerogel, there is obviously no 
size-dependence and we can safely conclude that the ad- 
sorption is gradual in the thermodynamic limit. On the 
contrary, in the lighter aerogel, the isotherms look steeper 
and steeper as the system size is increased and the re- 
sults suggest that there is indeed a discontinuity when 



L — > oo (the isotherm corresponding to L — 35 is shown 
here for completeness, but this size is probably too small 
to describe properly a 95% aerogel). As we discussed 
elsewhere [l^ in the case of fluid equilibrium behavior 
in purely random solids, one expects that the maximal 
slope of the isotherms should scale as L 3 / 2 at a first-order 
transition (assuming that the location of the largest jump 
fluctuates around its mean value ~p t (L) with a variance 
8/j, t (L) 2 oc L~ 3 ). As shown in Fig. 11, one can obtain 
a reasonably good collapse of the L = 50, 70 and 100 ir- 
reversible "compressiblity" curves dpf /dfi by using the 
scaling variable L x {p — ~p t (L)} with \ w 1.22. 

There are at least two possible explanations for the dis- 
crepancy with the expected value \ — 3/2. Firstly, the 
system sizes could still be too small so that too many gel 
realizations would have a non-typical behavior, with sev- 
eral avalanches of similar heights (note that in a coarse- 
grained picture of a 95% gel-fluid system where a re- 
gion of size £g is represented by a single effective spin, 
there would be only about 10 3 such effective spins in a 
sample of size L = 100). Secondly, one may be close 
to the critical value p c of the disorder, i.e., the critical 
value of the porosity for y = 2 and T/T c = 0.45, for 
which the infinite avalanche first appears. Studies of the 
T = RFIM have moreover shown that the critical re- 
gion is unusually large |T(|. The value x — 1-22 is com- 
patible with the predictions of Sethna and co-workers, 
with x — fi8/v ~ 2 — 77[3!|, although one cannot also ex- 
clude that the exponents differ from those of the conven- 
tional RFIM because of the presence of impurities j3||. In 
any case, it seems reasonable, on the basis of the present 
calculations, to predict the existence of a phase transi- 
tion with most probably a discontinuity in the adsorption 
isotherm in the thermodynamic limit, as shown in Fig. 
5(b). 
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FIG. 12: Desorption isotherms in 87% (a) and 95% (b) aero- 
gels. System sizes are L = 50 (a) and L = 100 (b). 



B. Desorption 

As discussed elsewhere [THj. fluid desorption in disor- 
dered porous solids can take place via several differ- 
ent mechanisms, depending on the temperature and on 
the structural and energetic properties of the solid (in 
Ref. ^| , however, only the influence of the interaction pa- 
rameter y, i.e., of the wetting properties of the adsorbed 
fluid, was described). In a first mechanism, desorption 
is due to the appearance of vapor bubbles in the bulk of 
the material, bubbles that grow, coalesce, and eventually 
extend over the whole pore space. The mass adsorbed 
then decreases continuously. In the other mechanisms, 
draining of the solid starts from the surface, and the 
desorption is associated to the penetration of a vapor- 
liquid interface which was previously pinned by the ir- 
regularities of the solid structure for p, larger than some 
threshold value /i c (the "depinning" threshold). The des- 
orption curve is then either gradual or discontinuous, de- 
pending on whether the growth of the vapor domain is 
isotropic and percolation-like (self-similar), or compact 
(self-affine). In those cases the transition is however al- 
ways critical because the scale of the rearrangements of 
the interface diverges as /i — > This is very simi- 

lar to the physics of fluid invasion in porous media[3^|. 
although we are considering here a single compressible 
fluid, and of field-driven domain wall motion in disor- 
dered magnets 38] . 

In order to find what are the relevant mechanisms in 
the 87% and 95% aerogels for y = 2 and T/T c = 0.45, 
we have studied the two systems in the presence and in 
the absence of the interface with the gas reservoir. In 
the latter case, we find that emergence of vapor bubbles 
in the bulk of the solid only occurs at low values of the 
chemical potential (fi ~ —5.27 and fi ~ —5.25 for 87% 
and 95% aerogels, respectively), very close to the value 
fJ'spi = —5.249 corresponding to the liquid (mean- field) 
spinodal of the bulk fluid at T/T c = 0.45. This shows 




FIG. 13: Invading vapor domain Vt (in dark) in a 87% aerogel 
at A* = -4.706 (L = 100). 10% of the fluid has drained out. 
The gas reservoir located at the bottom of the box and the 
aerogel are not shown. 



that the perturbation induced by the solid is too small 
to displace significantly the liquid spinodal. On the other 
hand, the contact with the ambiant vapor at the surface 
of the gel has a major influence, as illustrated in Fig. 12 
by typical desorption isotherms calculated in the presence 
of a free surface. One can see that when decreasing the 
chemical potential all the curves exhibit a pronounced 
drop much before vapor bubbles appear in the bulk of the 
material. This is a clear indication that the mechanism 
of desorption is due to the surface. 

The isotherms shown in Fig. 12 consist of many steps 
of varying size, but it also appears that one step is sig- 
nificantly larger than the other ones in most of the 95% 
samples, a feature that is not present in the isotherms of 
the 87% aerogel. This suggests that the growth morphol- 




FIG. 14: Same as Fig. 13 but for a 95% aerogel at fj, = -4.220. 
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ogy of the invading vapor domain may change with the 
porosity. This is confirmed by the snapshots displayed 
in Figs. 13 and 14 that show the invading vapor domain 
in two typical 87% and 95% samples when about 10% 
of the fluid has drained out of the aerogel. This is just 
after the onset of the sharp drop in the isotherms. As in 
the case of adsorption, we find that the distribution of 
the p^s is bimodal in this range of chemical potentials, 
with most of the fluid sites having a density lower than 
0.05 or larger than 0.95. One can thus identify unam- 
biguously the emptied sites. Although the volume occu- 
pied by the vapor is the same in the two samples (about 
200000 sites), the morphology of the domain is remark- 
ably different. In the 87% aerogel, the vapor domain 
exhibits an intricate isotropic structure that resembles 
that produced in invasion percolation. In contrast, the 
domain looks compact with a self-affine interface in the 
lighter aerogel. This is very similar to the two regimes 
observed in the T = RFIM when an interface sepa- 
rating two magnetic domains is driven by an external 
field [4(|: when the disorder is large, the interface forms 
a self-similar pattern with a large-scale structure char- 
acteristic of percolation, whereas the growth is compact 
and the domain wall forms a self-affine fractal surface at 
intermediate degrees of disorder (note that the use of the 
body-centered cubic lattice has allowed us to suppress the 
faceted growth regime that is observed at T/T c = 0.45 in 
the 95% porosity aerogel on the simple cubic lattice[3^|; 
this regime is an artifact of the lattice description) . 

In order to confirm the existence of two growth regimes 
and to determine the actual behavior in the thermody- 
namic limit, a finite-size scaling analysis of the desorption 
isotherms is required. Average isotherms calculated for 
different system sizes are shown in Fig. 15. By analogy 
with the problem of a driven interface in the T = RFIM, 
one expects that the total volume Vt(/Lt) of the invading 
vapor domain shows a power law divergence at the depin- 
ning threshold fi c . Then, assuming that the only relevant 
length scales near p, c are the sytem size L and a single 
correlation length £ ~ [([i — ^ c )/^ c ] _ly , the dependence 
of Vt on sytem size should be described by the scaling 
form L D f g(x) where Df is the fractal dimension charac- 
terizing the domain and g is a universal function of the 
scaling variable x — L x l v {yi — u„)/u„[40|. 

In both aerogels, one observes important finite-size ef- 
fects (see Fig. 15), but, unfortunately, they are not only 
due to the existence of a diverging correlation length in 
the system but also to boundary effects. There is indeed 
an initial regime in which desorption is due to the drain- 
ing of large crevices at the surface of the gel where the 
fluid is in direct contact with the ambiant vapor. For a 
given porosity, the number of these crevices is propor- 
tional to L 2 and the contribution to the fluid density is 
thus proportional to l/L. It turns out that this initial 
regime extends to rather low values of the chemical po- 
tential (p ~ -4.6 and fi ~ -4.15 for the 87% and 95% 
aerogels, respectively), so that the scaling region around 
the depinning threshold is too small to be studied prop- 
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FIG. 15: Average desorption isotherms in 87% (a) and 95% 
(b) aerogels for different system sizes. The number of gel 
realizations ranges from 1000 (L = 25) to 100 (L = 100) for 
the 87% aerogel, and from 1000 (L = 35) to 200 (L = 100) 
for the 95% aerogel. The arrow in (b) indicates the common 
intersection of the curves. 



erly and to extract the values of the critical exponents 
(such an effect is not present in the numerical studies of 
domain growth in the standard T — RFIM[3£| because 
there is a different random field on each site of the lat- 
tice and no equivalent of the crevices). We note, however, 
that the curves in Fig. 15(b) have a common intersection 
at /j, ~ —4.24, in contrast with those in Fig. 15(a). This 
is consistent with the scaling ansatz for the volume of the 
invading vapor domain with Df — 3 for the 95% aerogel 
and Df < 3 for the 87% aerogel. We thus conclude that 
the desorption seems to be discontinuous in the first case 
and gradual (percolation-like) in the second one. 

This leads to the isotherms shown in Fig. 5. Their 
shape resembles that of the experimental curves in Figs. 
4(b) and 4(c) of Ref. 0. The minor differences can be ra- 
tionalized: the experimental isotherm in the 87% aerogel 
(Fig. 4(c) in Ref.0) does not exhibit a sharp kink at fj, c , 
but this rounding may be due to the activated processes 
that are neglected in the present treatment; Fig. 4(b) in 
Ref. 8] corresponds to a 98% aerogel, which probably ex- 
plains why the shape of the hysteresis is more rectangular 
than in the present Fig. 5(b). 



C. Equilibrium 

In order to determine the equilibrium isotherms, one 
has to find, for each value of fx, the lowest lying state(s) 
among all the metastable states obtained from Eqs. 
(3). A complete enumeration of these states is, how- 
ever, an impossible numerical task, and like in previous 
work p^. Il4j |. we have only calculated a limited num- 
ber of states that, hopefully, can provide a good ap- 
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FIG. 16: Typical desorption scanning curves (points) and 
equilibrium isotherms (solid lines) in 87% (a) and 95% (b) 
aerogels (L=70). 




FIG. 17: Check of thermodynamic consistency along the ad- 
sorption, desorption, and equilibrium isotherms for a 87% 
porosity aerogel (L = 70). (a) adsorption and desorption 
(b) equilibrium. Symbols: average fluid density p/. Dashed 
curves: related quantity obtained by differentiating the cor- 
responding grand potential with respect to the chemical po- 
tential. 



proximation of the equilibrium isotherm. This can be 
checked a posteriori by using the Gibbs adsorption equa- 
tion Pf = —d(Cl /N) /dn which is only satisfied by the 
equilibrium curve [Taj. In Refs.fl^. Ibif. we searched for 
metastable states inside the hysteresis loop by starting 
the iterative procedure with initial configurations corre- 
sponding to uniform fillings of the lattice (with pf^ = p 
varying between and 1— p). This method, however, does 
not converge in dilute aerogels because the structure is 
very inhomogeneous and correlated, at least for distances 
smaller than £g . It is then likely that all metastable fluid 
configurations are also very inhomogeneous and cannot 
be reached iteratively from initial uniform fillings. We 
have thus searched for metastable states by calculating 
desorption and adsorption scanning curves, i.e., by per- 
forming incomplete filling or draining of the aerogel and 
then reversing the sign of the evolution of the chemical 
potential. For convenience, these (very long) calculations 
are done in systems with periodic boundary conditions in 
all directions, i.e., in the absence of an interface with the 
reservoir |4lj. 

Some typical desorption scanning curves are shown in 
Figs. 16(a) and (b). In both cases, the top curve is 
the major desorption branch, and by comparing with the 
isotherms shown in Figs. 12 or 15, it is clear that the 
part of this branch extending to low p's and correspond- 
ing to liquid-like metastable states that are isolated from 
the other states is an artifact coming from the absence 
of the interface with the reservoir. We also notice that 
there are no scanning curves in the upper part of the hys- 
teresis loop for the 95% sample. This is because there 
is a jump in the adsorption isotherm and, thus, there 
are no intermediate states from which desorption can be 
started. Therefore the states contained in this region of 
the plane (pf,p) cannot be reached by performing scan- 
ning trajectories[42j. (Further work is clearly needed to 



determine if there are no metastable states at all in this 
region or if the states can be obtained by other means, 
for instance by changing the temperature.) 

Fig. 16 also shows the approximate equilibrium 
isotherms obtained by selecting, for each value of p, the 
state a that gives the lowest value of the grand potential, 
as calculated from Eq. (4). We have checked that taking 
into account the additional metastable states obtained 
from the adsorption scanning curves did not change sig- 
nificantly the results (as in Ref. , we have also checked 
that keeping all solutions with a weighting factor equal to 
the Boltzmann factor gives the same isotherms). As illus- 
trated in Figs. 17 and 18, the Gibbs adsorption equation 
is very well verified along these curves. This indicates 
that we have indeed obtained a good approximation of 
the true equilibrium isotherms. In contrast, thermody- 




(b) 
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FIG. 18: Same as Fig. 17 but for a 95% aerogel. 
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FIG. 19: Average equilibrium isotherms in 87% (a) and 95% 
(b) aerogels for different system sizes. The number of gel 
realizations ranges from 500 (L = 25) to 20 (L = 100) for the 
87% aerogel, and from 2000 (L = 35) to 200 (L = 100) for 
the 95% aerogel. 



namic consistency is strongly violated along the adsorp- 
tion and desorption branches. The presence of (delta) 
peaks in dVl/dp also shows that the grand potential f2 
changes discontinuously (as the fluid density pf) during 
adsorption and desorption in a single finite sample. On 
the other hand, Q is continuous (but pj is discontinu- 
ous) al ong the equilibrium isotherm, as already noticed 
in Ref.[l|. 

One can see from Figs. 16(a) and 16(b) that the equi- 
librium behavior is different in the 87% and 95% poros- 
ity aerogels. In particular, there is a large final jump 
in the isotherm of Fig. 16(b). The same feature exists 
in all samples, but in order to conclude on the actual 
behavior in the infinite system one has again to analyze 
the size-dependence of the average isotherms. This is 
shown in Fig. 19. For the 87% aerogel, there is al- 
most no size-dependence and it is clear that no transi- 
tion occurs when L — > oo. On the other hand, for the 
95% aerogel, the isotherms become steeper as L increases 
and there is a rather well defined common intersection at 
/i t ~ —4.17. This is a clear indication that a phase tran- 
sition does occur in the thermodynamic limit. However, 
we have not succeeded in obtaining a satisfactory col- 
lapse of the isotherms using the scaling reduced variable 
L 3 / 2 [/x — nt(L)]/ nt(L) as was done in Ref.|14j in the case 
of a random solid. The origin of this problem is still 
unclear, but we suspect that the system sizes are again 
too small. Indeed, we note that the maximal slope in- 
creases rather weakly between L — 35 and L = 70, but 
has a size-dependence consistent with the exponent 3/2 
between L = 70 and L = 100. We thus conclude that 
our results are consistent with a discontinuous jump in 
the thermodynamic limit, as indicated in Fig. 5(b). This 
corresponds to a true equilibrium liquid- vapor phase sep- 
aration in the system. 



V. SUMMARY AND CONCLUSION 

In this paper, we have proposed an interpretation of 
the hysteretic behavior observed in the experiments of 
4 He adsorption in light silica aerogels. The overall shape 
of the experimental hysteresis loops is well described by 
our theoretical model and we have been able to reproduce 
the dramatic influence of porosity. In this interpretation, 
the disordered character of the aerogel structure plays an 
essential role, the porosity p being the tunable parame- 
ter that controls the amount of disorder in the system. 
The history-dependent behavior is thus associated to the 
presence of many metastable states in which the system 
may be trapped and that prevent thermal equilibration at 
low temperatures. The most important conclusion of our 
study is that adsorption and desorption obey to different 
mechanisms and may be gradual or discontinuous, de- 
pending on the porosity. Adsorption is insensitive to the 
presence of the interface between the solid and the exter- 
nal vapor, and the change in the shape of the isotherm 
from smooth to discontinuous as p increases from 87% 
to 95% has been related to the appearance of a infinite 
avalanche occuring in the bulk of the system, similar to 
what happens in the T = RFIM below the critical 
disorder 1 16j. In contrast, desorption is triggered by the 
presence of the outer surface and an associated depin- 
ning transition. The change in the shape of the isotherm 
is then related to a change in the morphology of the in- 
vading vapor domain from percolation-like to compact 
(note however that the mechanism may be different in 
aerogels of lower porositv|43j). 

One must keep in mind that our calculations are based 
on local-mean field theory in which disorder-induced fluc- 
tuations are properly accounted for but thermal fluctua- 
tions are neglected. This appears to be a reasonable as- 
sumption at very low temperatures where the time scale 
associated to thermally activated processes is much larger 
than the experimental time scale. In this case, both 
the adsorption and desorption branches are metastable 
and we have shown that the true equilibrium isotherm 
(which probably cannot be reached experimentally) is 
somewhere in between: this isotherm also changes from 
gradual to discontinuous as p increases. The situation 
is different at higher temperature (in the vicinity of T c ) 
since true equilibrium behavior without hysteresis has 
been observed experimentally 0, 0- It will be therefore 
of interest to study the influence of temperature on the 
hysteretic behavior and to understand how the dynamics 
of relaxation towards equilibrium changes with T (see, 
e.g., Ref.^3] f° r a recent study of the relaxation behav- 
ior associated to capillary condensation in a lattice model 
of Vycor) . 

Finally, we suggest that the scenario for filling and 
draining described in this work could be tested in more 
detail by performing experiments in a series of aero- 
gels of gradually increased porosity. If our interpreta- 
tion of the adsorption process in terms of avalanches is 
correct, there should exist a value of the porosity for 
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which the isotherm becomes critical and the distribu- 
tion of avalanche sizes follows a power law with well- 
defined critical exponents|l6j|. Using superfluid instead 
of normal 4 He could perhaps allow to observe distinct 
avalanche events in the aerogel and to study their sta- 
tistical properties. Such a study has been performed re- 
cently in the nanoporous material Nucleopore by a ca- 
pacitance technique 45] . One could also check that the 
draining of the aerogel starts from the surface and that 
the growth of the invading vapor domain obeys different 
regimes. Related studies have been for instance carried 
out in Vycor using ultrasonic attenuation and scattering 
techniques |32l |46|. Such experiments would give a defi- 



nite answer to the long-standing question about the na- 
ture of hysteresis in fluid adsorption in disordered porous 
media. 
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